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We study the quark mass dependence of the QCD phase transition by an effective potential defined 
through the distribution function of observables. As a test of the method, we study the first 
order deconfinement phase transition in the heavy quark mass limit and its fate at lighter quark 
masses. We confirm that the distribution function for the plaquette has two peaks indicating that 
the phase transition is of first order in the heavy quark limit. We then study the quark mass 
dependence of the distribution function by a reweighting method combined with the hopping 
parameter expansion. We find that the first order transition turns into a crossover as the quark 
mass decreases. We determine the critical point for the cases of Nf = 1, 2, 3 and 2 + 1. We find 
that the probability distribution function provides us with a powerful tool to study the order of 
transitions. 
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1. Introduction 

It is of great importance to know the nature of the QCD phase transition in understanding the 
evolution of the early universe. The nature of the phase transition depends on the quark masses. 
The deconfinement phase transition is expected to be first order when the masses of up, down and 
strange quarks are either sufficiently large or small, and is crossover in the intermediate region be- 
tween them. At the physical quark masses, previous studies with staggered quarks strongly suggest 
that the transition is crossover [[TJ, ||]. A confirmation of this result by other fermion formulations 
such as Wilson-type fermions or by other methods for the analysis, however, is mandatory to draw 
a definite conclusion on the order of the transition at physical quark masses. 

In this report, we propose a method based on measurements of the probability distribution 
function of a physical quantity to identify the order of the phase transition. Since two phases 
coexist at a first order phase transition point, one can identify the first order phase transition by 
measuring the distribution function, which must have two peaks at the first order phase transition 
point where two phases coexist with equal probability. In our study we take the plaquette, i.e. lxl 
Wilson loop, to label the state. 

As a test of the distribution function method, we investigate the quark mass dependence of the 
order of the QCD phase transition in a large quark mass region, where numerical calculations are 
much easier than those in the light quark mass region. In addition, the reweighting [Q, |j] technique 
is used to vary quark masses. 

This paper is organized as follows: Basic properties of the plaquette distribution function are 
discussed in Sec. ||[ The method to calculate the plaquette distribution function by the hopping 
parameter expansion in the heavy quark mass region is introduced in Sec. |3[ In Sec. [|, we present 
results for an effective potential defined from the distribution function, and show that the order of 
phase transition changes from the first order to crossover as the quark mass decreases from infinity. 
We then evaluate the location of the critical point. The paper is summarized in Sec. H. 



2. Probability distribution function 



The probability distribution function provides us with one of the most fundamental approaches 
to identify the order of the phase transition. Because there exist two phases simultaneously at a first 
order phase transition point, we expect that the probability distribution has two peaks there. In this 
report, we study the distribution function of the average plaquette P, i.e. lxl Wilson loop. We 



use the plaquette action (S g ) 
quark part: 



for the gauge part and the standard Wilson quark action (S q ) for the 
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where = N% x N t is the number of sites, Nf is the number of flavors, Kf is the hopping pa- 
rameter, and j3 = 6/g 2 . The quark mass is controlled by Kf, which is proportional to when 
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Kf is small, while the lattice spacing is mainly controlled by /3. For the case of degenerate quark 
masses, i.e. Ky = K for / = 1, • • • ,Nf, the probability distribution function for the plaquette variable 
is defined by 

w{P',f3,K) = J ®U 8{P'-P[U]) (detM(K)) Ni e 6 P N ^ p , (2.4) 

where 8(x) is the delta function. The partition function is given by 3f(K,fi) = J w(P',/3, K)dP' . In 
the folio wings, we denote P' simply as P. 

The plaquette distribution function can be obtained by the histogram of P and has the following 
useful property: Under the parameter change from /$o to /3, the weight w(P,[3,k) becomes 

w{P,p,K) = e 6 ^- MN ^ p w(P,po, K). (2.5) 

The effective potential defined by 

V(P,P,K) = -\nw(P,P,K) (2.6) 

varies as 

V(P,P,K) = V{P,fa, K) - 6(0 - j3o)AT site P (2.7) 
under this parameter change. From this property, we find that dV / dP changes trivially as 

dV dV 

— (P,P, K) = ^ (P, A), k) ~ 6(J8 - ftRite, (2.8) 

and d 2 V / dP 2 is unchanged under the j3 shift. Therefore, the shape of dV / dP as a function of P 
does not change with /3 up to the j3 -dependent constant. 

When a first order transition exists, the distribution is a double-peaked function at the transition 
point. The effective potential is then a double- well function and the derivative of V is an S -shaped 
function. In the transition region, dV jdP vanishes at three points. To locate the transition point, 
the fine turning of j8 to the transition point is not required because the P-dependence of dV / dP 
itself is independent of /3 . Therefore, the measurement of dV jdP is useful to identify the first order 
phase transition. 

Next, we discuss the fc-dependence of V considering the ratio of the distribution functions at 
different K and Kq\ 

( 8(P'-P[U]) W*)f* \ 

pfpi „ „ n = w(P',p,K) _ f®U8(P'-P{U])(d e tM(U,K)) N f _ \ (dctM(t/.«a)) f/p^ 

' - MP' ~ f®U8(P>-Piumd e tM(U,Ko)) N f ~ (8(P'-P[U])) {fStKo) ■ ^> 

This R(P, k, Kb) is independent of j3. Using R(P, K, Kb), the K-dependence of the effective potential 
is given by the following equation: 

V(P,P,K) = - \nR(P, K, Kb) + V(P,P,Ko), (2.10) 

under the change from Kb to K. 

This argument can be generalized to improved gauge actions easily. For the case of im- 
proved gauge actions including larger Wilson loops, we should define the average plaquette as 
P = —Sg/ (6Af s itej3), which is a linear combination of Wilson loops. On the other hand, /3 -dependent 
improved quark actions make the analysis more complicated. 
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3. Quark determinant in the heavy quark mass region 



To investigate the quark mass dependence of the plaquette effective potential, we evaluate the 
quark determinant by the Taylor expansion with respect to the hopping parameter K in the vicinity 
of the simulation point fc : 



In 
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detM(Kb) 
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Calculating the derivative of the quark determinant S> n , the ^-dependence of the effective potential 
can be estimated. 

In this study, we focus on the boundary which separates the first order transition region near 
the quenched limit and the crossover region. The boundary is expected to exist near k = 0. We 
adopt K"o = where M xy = 8 X} y. The (dM / 'dK) xy is the gauge connection between x and y, and the 
nonzero contribution of & n is given by Wilson loops and Polyakov loops. Considering the anti- 
periodic boundary condition and gamma matrices in the hopping terms, the leading contributions 
of the Taylor expansion are given by the K 4 term and k n ' term: 



In 



detMf k) 



detM(O) 
where Q. is the Polyakov loop: 



28 8AW K A P + 1 2 x 2 N 'N^ fc^' Reft + • • ■ 
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^n,4^ n+ 4 j4 ^ n+ 24,4 • • • U n +{N t - 1)4,4 



(3.3) 



(3.4) 



The ratio R(P, K,0) = w(P,j8, Jc)/w(P,j8,0) is calculated by the following equation in the region of 
small K for any j3 : 

. l (d(P'-P[U})exp[NfN^(l2x2 N 'K N 'Rea[U} + ---)]) m ^ m 
R(P',K,0) = e N ^ 2 ^ p — 1 U „ rTrlxs — n/ V } > K0=O) .(3.5) 



(8(F-P[U])) 



(j8,Kb=0) 



For N t = 4 the truncation error is 0(k 6 ). The contribution from the plaquette can be absorbed in 
the redefinition of /3. Using Eqs. ( |2.10| ) and (3^5), we investigate the K"-dependence of the effective 
potential in the heavy quark mass region. 



4. Numerical simulations and the results 

In the heavy quark mass limit we perform simulations of SU(3) pure gauge theory on a 24 3 x 4 
lattice. To generate the configurations, the pseudo heat bath algorithm of SU(3) gauge theory is 
used, and the over relaxation is performed 4 times every update. Because the effective potential 
must be investigated in a wide range of the plaquette value to identify the order of the phase tran- 
sition, we performed simulations at five points in the range of j8 = 5.68 - 5.70. In the calculation, 
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Figure 1: Derivatives of the effective potential for 
pure gauge at /3 =5.68-5.70. The black line is the 
average. 
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Figure 2: The K"-dependence of the derivative of 
the effective potential at /3 = 5.69. 



we approximate the delta function by a Gaussian function: 8{x) m 1 /(A-^/n) exp [— (jc/A) 2 ]. Ex- 
amining the resolution and the statistical error in the distribution function, we adopt A = 0.000283. 

We then calculate the derivative dV/dP by the difference between the potentials at P and 
P + AP. We adopt AP = 0.0001 with which the AP-dependence is much smaller than the statistical 
error. The results of dV jdP at K = are shown in Fig. [jj In this figure, we adjust results at 
different j3's to j8 = 5.69 by using Eq. (2.8). The results of dV jdP obtained in simulations at 
different j8 is consistent within errors, though the ranges of P in which V is reliably obtained are 
different. The jackknife method is used to estimate the statistical error of the effective potential and 
its derivatives. We combine these data obtained by an weighted average with the inverse-square of 
each error, which is shown by a black line in Fig. [I]. The data is excluded from the average over 
different ensembles if the error is large and P is far away from the peak of the distribution function 
at each j3 . We find, from this figure, that dV jdP at K = is not a monotonically increasing function, 
meaning the double-well structure of the effective potential. 

The quark mass dependence of the effective potential is investigated by calculating R(P, K,0) 
at the order of K 4 in Eq. (j33[). Using the data of w(P,fi,K = 0) and R(P, K,0), we evaluate the 
dV jdP up to a finite value of fc. The results for Nf = 2 are plotted in Fig. ^[ The S-shape structure 
becomes milder as fc increases and seems to turn into a monotonically increasing function around 
fc = 0.066. This behavior suggests that the first order phase transition at K = becomes weaker as 
K increases and the transition changes to crossover at fc « 0.066 or larger. 

To evaluate more precisely the value of K at the boundary where the first order phase transition 
is terminated (fc cp ), we calculate the second derivative of V(P,fi, k) by a numerical differentiation 
of dV jdP. When the transition is first order, there is a region where the second derivative of 
V(P,j8, k) is negative between the two bottoms of V(P,P,k). The first order transition region is 
thus identified by measuring the sign of the second derivative of V. As discussed in Sec. ||, the 
second derivative of V(P,fi, k) is independent of j8. Therefore, the identification can be performed 
without considering the j3 -dependence. 

We plot the results of d 2 V jdP 1 at K = 0.058, 0.062 and 0.066 for N f = 2 in Fig. |. The range of 
P where d 2 V jdP 2 < becomes narrower at non-zero K than that at K = 0. We calculate d 2 V /dP 2 
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Figure 3: The second derivative of the effective potential at K = 0.058,0.062,0.066. The black symbols are 
the results at K = 0. 
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Figure 4: The k value at which d 2 V /dP 2 = 
for each fixing P for Nf = 2. 
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Figure 5: The quark mass dependence of the order of 
phase transition in 2 + 1 flavor QCD. 



varying fc with fixing P, and find the value of fc at which d 2 V/dP 2 = 0. The result is plotted in 
Fig. ^ for Nf = 2. In the region below the symbols, the curvature of V is negative. Fig. ^ shows a 
peak around P = 0.5478, and d 2 V /dP 2 is positive for all P when fc is larger than the peak height 
0.068(7). This means that the phase transition is no longer first order above the critical value of 
Kc P = 0.068(7) for Nf = 2. We calculate the maximum value of fc where d 2 V jdP 2 = , which is 
Kcp, with estimating the jackknife error. The error is estimated by a jackknife method. 

This analysis can be also applied to Nf = 2 + 1 QCD having different light quark mass and 
strange quark mass as well as Nf = 1 and 3 cases. For the case of degenerate quark masses, the 
second derivative of V does not change if the term NfK N ' is constant. From this property and the 
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In 



288AU(2f4j + 4)P + 12 x 2 N 'N*{2^ + )Reft + • • • ,(4. 1) 



result of the critical point for N f = 2, we find jq p = 0.081(8) for N f = 1 and 0.061(6) for N f = 3. 
The result of Nf = 1 is consistent with the results obtained by an effective Z(3) model in Ref. [g]. 
For the case of A/f = 2 + 1, within the leading order of the hopping parameter expansion, the quark 
determinant in the partition function is given by 

' (det M ( fc ud ) ) 2 det M ( fc s ] 
_ (detM(O)) 3 

where K U( \ and K s are hopping parameters for light and strange quarks. Because the contribution 
from the plaquette term in this equation does not affect the second derivative of V, the difference 
from the case of Nf is just the replacement from NfK N ' to 2fc^!j + fcf'. Thus the line which separates 
the first order phase transition and the crossover is given by 

2«2 + K* = (k3= 1 )*, ( 4 - 2 ) 
where N t = 4 and K^p 1 = 0.081(8) in this study. We draw the line in Fig. |5[ 

5. Conclusion 

We studied the order of the deconfinement phase transition in the heavy quark mass region of 
QCD by calculating the probability distribution function of the average plaquette. The distribution 
function at K = is evaluated by measuring the histogram in quenched QCD, and the ^-dependence 
is investigated by using the reweighting method within the approximation by the leading approxi- 
mation of the hopping parameter expansion for IndetM. Because the distribution function must be 
measured in a wide range of P to study the order of the phase transition, we performed quenched 
simulations at 5 points of j8 and combined these results. 

We found that the distribution function shows two peaks at K = indicating the first order 
transition, and the double-peaked shape becomes weaker as K increases. We calculated the critical 
point Kc p . The results are JQp = 0.081(8), 0.068(7) and 0.061(6) for N f =1,2 and 3, respectively. 
The boundary of the first order transition region for Nf = 2 + 1 QCD is also determined in the 
two-dimensional parameter space of K u a and fc s . 
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